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Abstract. A connection is made between the exact eigen states of the BCS Hamiltonian and 
the predictions made by the Tamm-Dancoff Approximation. This connection is made by means 
of a parametrised algebra, which gives the exact quasi-spin algebra in one limit of the parameter 
and the Heisenberg-Weyl algebra in the other. Using this algebra to construct the Bethe Ansatz 
solution of the BCS Hamiltonian, we obtain parametrised Richardson-Gaudin equations, leading 
to the secular equation of the Tamm-Dancoff Approximation in the bosonic limit. An example 
is discussed in depth. 



1. Introduction 

The description of superconductive properties in many-particle systems such as metals pQ or 
atomic nuclei [2j involves the process of pairing. Particles in time-reversed single-particle states 
will lower their energy by forming a pair, constituting a condensate in the ground state of the 
system. This ground state is very-well approximated by the Bardeen-Cooper-Schrieffer (BCS) 
wavefunction [lj, which is a variational solution of the BCS Hamiltonian in the grand canonical 
ensemble. While this solution is essentially sufficient for the thermodynamic limit, N — > oo, it 
is not fit for the description of finite-size systems. This is due to the increasing relevance of the 
pair fluctuations with decreasing number of particles, which cannot be taken into account in the 
BCS approximation [3]. This is not of paramount importance for modestly small systems, such 
as e.g. single-species nucleon pairing within a single major nuclear shell [3], because these can 
be diagonalised numerically on a modern computer with relative ease. However, for mesoscopic 
systems, such as metallic nanograms [5], or nucleons in multiple open shells, the dimension of the 
Hilbert space is huge, such that straightforward numerical diagonalisation becomes practically 
intractable. 

Fortunately, alternative methods to extract the spectroscopy of the model (BCS) Hamiltonian 
exist. Richardson and Sherman [BJ [7] have shown that the reduced BCS Hamiltonian, with a 
level-independent interaction, is exactly solvable by means of a Bethe Ansatz product wave 
function, provided the variables in the Ansatz wavefunction provide a solution of the set of non- 
linear Richardson-Gaudin (RG) equations. As a result, the problem of finding the eigenvalues 
and eigenstates of the BCS Hamiltonian is reduced from diagonalising a matrix that scales 
factorial with the size of the system to solving a set of algebraic equations that scales linear 



with the system size. However, the RG equations are highly non-linear in the variables, and 
proved to be very challenging to solve 0(9], due to the occurrence of singularities at the critical 
interaction values. It was only recently that this obstacle has been sufficiently removed [TO] 
by means of a change of variables method, canceling the singularities at the critical points. 
Alongside the numerical usefulness, the method also contributed to a better understanding of 
the critical points [11] , which finds its place among other recent studies of the structure of the 
eigenstates of the exactly solvable and integrable models [1'2\ [TO] . For an extensive review on 
the subject, I would like to refer to the recent colloquium [14| . 

The purpose of the present manuscript is to contribute to the understanding of the Bethe 
Ansatz eigenstates of the BCS Hamiltonian and the associated RG equations by making a 
connection with the Tamm-Dancoff Approximation (TDA) [TO]. In the next section, the BCS 
model and its algebraic properties will be introduced. In section O the TDA is introduced and 
its relation with RG is made. Conclusions are presented in section [J] 

2. The reduced BCS model, an exactly solvable model 

Consider a fermionic system with n single-particle energy levels (1 < i < n) with degeneracy 
f2j. For condensed matter physics the levels are two-fold degenerate (f2j = 2), corresponding 
to the two projections of the electron spin. In the context of nuclear structure physics, we can 
assign a half-integer total angular momentum ji with degeneracy f2 i = 2ji + 1 to every level [1] . 
Naturally, and in order to make the discussion as general as possible, we will treat the case of 
condensed matter as ji = ^ systems. The reduced BCS Hamiltonian for such a general system 
is given by 

n n 

H = ^2ein i + g^2 SjS k , (1) 

i=l i,k=l 

with the operators in the Hamiltonian given by 

k k 

mi=-ji rrii=-ji 

which are known respectively as the particle counting, pair creation and pair annihilation 
operator. The underlying algebraic structure of these operators is the su{2) quasi-spin algebra. 
This algebra is spanned by the pair creation/annihilation operators, completed with the operator 
Sf = \hi - lOj. We obtain 

[S?,sl] = 5 ik s}, lS°,S k ] = -5 ik Si, [sj,S k ]=25 ik Sf. (3) 

So, the spectrum generating algebra of the pairing problem is ©™ =1 su(2)i, with su(2)i the quasi- 
spin algebra for the level i. The associated Hilbert space is constructed from the direct product 
of the physically allowed su{2)i representations m s .), with the quantum number^ defined 
by the eigenvalues of the quadratic Casimir operator Cf of SU{2)i and its SO(2)i projection Sf 

C?\si,m Si ) = Si(si + i)\si,m Si ) (4) 
Sf\si,m Si ) = m Si \si,m Si ). (5) 

The physically allowed representations are those for which the following relations hold 

Ttl Si — 2^"* 4^i> &i — 2^*' 



1 I will denote the quasi-spin projection quantum number m of the quasi-spin Si explicitly by m Bi to rule out 
any confusion with the total angular momentum projection rrii of ji. 



with m the number of particles in a level i and Vi the seniority of the level i, which is the number 
of particles in a level i not pairwise coupled |16| . Obviously, for systems with a definite total 
number of particles N, we require the extra condition Y2% n % = N. For simplicity, we will assume 
that there are no unpaired particles present in the system (wj = 0, Vi), so the total number of 
pairs is given by N/2. 

The key ingredient of the Richardson-Gaudin treatment of the pairing model is that the 
Hamiltonian ([T]) can be diagonalised by means of a Bethe Ansatz wavefunction [H] 

N/2 / n „ t \ 
o=l \i=l / 

provided the Richardson variables, given as the free parameters E a , are a solution of the RG 
equations 

n N/2 

1 + 2 ^E^V-^E^^ = ' (l<«<iV/2). (8) 

1=1 fij^OL r 

The energy of the eigen state is then given by E = Yla=i A*- In the next section, it 
will be discussed how the Richardson-Gaudin treatment can be related to the Tamm-Dancoff 
Approximation, which is a well-studied method in many-body physics. 



3. Richardson-Gaudin and the Tamm-Dancoff Approximation 

3. 1 . The Tamm-Dancoff Approximation 

Most commonly, the Tamm-Dancoff Approximation is used in the context of particle-hole 
excitations across a closed Hartree-Fock vacuum |15[ S] , but we can interprete it here somewhat 
larger as the description of elementary excitations which are also related to the Random Phase 
Approximation or Equations of Motion methods [T7]. The basic idea behind the approximation 
is that we can treat all the excited states of the Hamiltonian as harmonic excitations or 
multiphonon states of the elementary eigen modes or phonons of the system. In the case of 
the BCS Hamiltonian (JT]), the elementary excitations can be found by solving the eigen value 
equation for the 1-pair excitations b\ = Y^i YkiSj 

Hbl\9) = huj k bl\e). (9) 
Once this is solved, the excited states are approximated by the following multiphonon states 

n 

fc=i 

with = N/2 and the notation [v^\ is short-hand for the array [z^i, 1*2, ... , Vn\- Similarly, 

the energy of this state is equal to 

n 
k=l 

Obviously, this is a gross simplification compared to the exact Bethe Ansatz wave function 
(J7]). However, it is worthwhile to scrutinise the differences. For a 1-pair system, the RG and 
TDA approach are identical because the action of the elementary phonons b\ on the vacuum is 
equivalent to the Bethe Ansatz wavefunction (|7|) for a single pair {N/2 = 1). Moreover, it is 



Table 1. The parameters of a simple 3-level pairing system that serves as an example throughout 
the text. 



a 



Si 



I 6 0.0000 
J 4 0.7796 
| 10 2.1024 



well-known that the TDA eigen modes (huik) for the reduced BCS Hamiltonian can be identified 
as the solutions of the secular TDA equation [H [J5] . 

i=i 3i 

which is identical to the RG equation for a single pair ([8]). Naturally, the situation is different 
for multi-pair excitations. It is instructive to continue with an example. Consider a system of 
4 pairs in 3 levels with the parameters given in Table [TJ The left-hand side of the secular TDA 
equation (fl"2"|) for the parameter set of Table [Q is plotted in Figure [lj revealing quite a bit of 
information at a glance. It is readily understood that there are as many eigen modes as there 
are levels and one eigen mode is structurally different from the others in the sense that it is not 
bound in a domain between two single-particle poles. This solution is unbound below and is 
generally referred to as the collective solution. 

8 

6 

4 

2 


-2 
-4 
-6 
-8 

-6 -4 -2 2 4 6 
ho) [a.u.] 

Figure 1. The TDA secular equation (eq. (|12p ) for the parameters given in Table [1] and an 
interaction strength g = —0.5. The poles at the single-particle energies are accentuated by 
vertical lines and the 3 solutions of the equation are marked graphically by open circles. 




So, we gather that there are 3 eigen modes for the example under discussion. Since we are 
accommodating 4 pairs in the system, we can construct a total of 15 distinct TDA multiphonon 
states [z^i, U2, ^3] with i>i + V2 + v z = 4. Comparing this with the real dimension of the Hilbert 
space, which is 11, we immediately realise that TDA overestimates the total number of physically 
allowed states. The question arises which ones are physical and which are not. To investigate 



this question, it is instructive to plot the exact excitation spectrum of the system against the 
predictions made by TDA, which is done in Figure [2j We can roughly divide the Figure in 2 
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Figure 2. The exact spectrum of the BCS Hamiltonian (JTJ) in (black) solid lines and the 
predictions of the TDA in (blue) dashed lines are plotted for a range of interaction strength g. 
In the inset, we find a zoom of the plot for small values of g. 

different domains, i.e. the small-interaction domain \g\ ~ and the large-interaction domain 
\g\ S> 0. The former domain is enlarged and depicted in the inset of the Figure. From this inset, 
it is clear which TDA states are unphysical for the small-interaction domain. We see that for 
\g\ — > 0, the exact ground-state energy is given by approximately 6e\ + 2e2 = 1.5592, whereas 
the TDA predicts the ground state at 8s\ = 0.0000. This is because the TDA can put all 4 
elementary excitations in the lowest level, [^1,^2^2] = [4,0,0], whereas this is fundamentally 
forbidden in the Hilbert space because of the Pauli principle. The degeneracy of the lowest level 
Oi = 6 only allows for 3 pairs, so one pair needs to be promoted to the 2nd level. As a result, 
we can identify [4, 0, 0] as physically forbidden in this domain. Further investigation reveals that 
[1,3,0], [0,4,0] and [0,3,1] are also forbidden on similar grounds. The situation is less clear for 
the large-interaction domain because we cannot rely on pure combinatorial arguments. However, 
we can draw qualitative conclusions by noticing that the excitation spectrum tends to reorganise 
in small clusters around \g\ = 1. We count clusters of 1, 2 and 3 states from the bottom up and 5 
more states that have no cluster characterisation. Inspecting the predictions by TDA, we notice a 
similar pattern where every cluster can be associated with a distinct number of collective phonons 
([4,0,0] for the ground state, [3,1,0] and [3,0,1] for the next cluster, etc.). Nevertheless, we 
also notice that TDA systematically underestimates the excitation energy of the corresponding 
states. Again, this is a manifestation of the Pauli exclusion principle, however in a different 
way. The correct clustering for the lowest states reveals that TDA gives a correct idea about 
the collective structure of the lowest states, but TDA allows Pauli- forbidden configurations to 
contribute to the energy of the eigen state. Anyhow, these qualitative arguments do not answer 
the question which TDA states are related to which exact eigenstates and which TDA states are 
not physical. In the next subsection, I will present a construction for this particular purpose. 
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3.2. An algebraic approach to go from Tamm-Dancoff to Richardson- Gaudin 
We may conclude from the previous subsection that the difference between TDA and the exact 
solution are due to the Pauli principle. In algebraic phrasing, the fundamental algebra of the 
exact pairing model is quasi-spin su(2) whereas TDA is characterised by the bosonic Heisenberg- 
Weyl algebra hw(l), given by [3]. 

[b k ,bj] = s kl i, [S fc ,i] = [S£,i] = o. (13) 

From a physical point of view, one can conclude that the fermion pairs Sj (eq. ([2}) have been 
replaced by genuine bosons b\ in the TDA. This substitution is rather abrupt, so it would be 
interesting to have a more adiabatic means to bosonise the fermion pairs. This can be realised 
by introducing the following parametrised algebra 

[Sf, Si) = 8 ik sj, [Sf, S k ] = -S ik Si, [Sj, S k ) = M£2S° + (£ - 1)^4), (14) 

with £ a real- valued parameter between and 1. On the one hand, we retain the quasi-spin 
algebra for £ = 1, and on the other hand, we get the Heisenberg-Weyl algebra in the £ = limit, 
if we scale the boson creation- and annihilation operators with a factor W2/Qi, and interprete the 
operator Sf as a boson counting operator. In the intermediate region, we can use the parameter 
£ as a continuous switch to transform the fermion pairs into bosons and the other way around. 
This construction is closely related to the contraction of u{2) to hw(l), as defined by Arecchi 
et. al. [18\ \19\ . The question is now whether this functionality is also transferred to the pairing 
system. Remarkably, if we use the parametrised algebra (|14p instead of the original quasi-spin 
algebra ([3]), we observe that the system remains exactly solvable by means of a Bethe Ansatz 
wavefunction, provided the free variables in the wavefunction are a solution of the following 
parametrised RG equations 

n io _f. N/2 
l^Ef^f "^E^ = 0, U<*<AV2). (15, 

It is readily verified that this set of equations reduces to the correct corresponding expressions 
in the limits for £. For £ = 1, we retain the set of RG equations (eq. (|HJ)), and for £ = 0, we 
obtain N/2 identical copies of a (seniority zero (vi = 0, Vi)) secular TDA equation. This means 
in the latter case, that every individual pair in the Bethe Ansatz wavefunction needs to satisfy 
the TDA secular equation, independent from the other pairs because the coupling term between 
the variables E a in the eqs. (|15p has disappeared at £ = 0. More details on the technical aspects 
of these results, as well as on the nature of the algebra (|14p will be given in a forthcoming paper. 

Once we have the ^-parametrised RG equations, we can investigate how the exact eigen states 
of the pairing problem are related to the solutions of the TDA. For this purpose, we can start 
from a known solution of the pairing problem [10] and see which TDA solution corresponds to 
this solution, or vice versa, we can start from a particular multiphonon solution of TDA and 
reconstruct the corresponding exact solution, if it exists. This can be done in both cases by 
adiabatically switching off/on the £ parameter in the eqs. (|15p and solving numerically along 
the trajectory. One particular result for the present example (with parameters from Table [1] 
and g = —0.5) can be found in Figure [3j We start from the [4,0,0] multiphonon TDA solution 
of the system (see also Figure [1] for the TDA eigen modes) and, by adiabatically switching on 
the parameter £, obtain the solution of the RG equations. In the Figure, the ^-trajectory of the 
solution in the complex plane is depicted. 

Finally, we can use this method to make a one-to-one connection between the exact eigen 
states and the predicted states by TDA, plotted in Figure [2j This is done extensively in Figure H] 



-3.5 -3 -2.5 -2 -1.5 -1 
Re(£J 



Figure 3. The ^-trajectory of a solution of the parametrised RG equations (|15p in the complex 
plane. We start from the [4, 0, 0] TDA solution and adiabatically go to the corresponding solution 
of the RG equations ([8]) for £ = 1 (denoted by E a , a = 1 ... 4). The parameters for the model 
are given in Tableland g = —0.5. 



andO Each row of panels corresponds to a particular TDA state, denoted by [v\, zv 2 , z/ 3 ]. In the 
first panel, the energy of the elementary TDA eigen modes (blue dashed line) and the real part 
of the corresponding RG variables (black full line) are plotted as a function of the interaction 
strength \g\. The second panel shows the imaginary part of the RG variables, and the third 
panel shows the total energy of the exact eigenstate Y2 a =i E a (black full line) as well as the 
TDA prediction of the total energy (eq. (jlip) (blue dashed line) as a function of \g\. For every 
TDA state [i>i, z^, ^3], the interaction strengths g is scanned and the corresponding RG variables 
are given if they exist. In some cases, such as e.g. for the low-interaction domain of [4, 0, 0] (first 
row of Figured]), the TDA state did not produce an exact solution and it is only for values larger 
than \g\ ~ 0.190 that it leads to an exact state (the ground state). It is readily seen from the 
second row that the ground state below this value \g\ ~ 0.190 is connected to the [3, 1,0] TDA 
state. Dotted lines are added to the figures to guide the eye at those transitional points where 
the TDA origin of the exact states changes from one state to another. In addition, a label \i) 
is added in the third panel of every row to denote the corresponding exact state. We can draw 
multiple conclusions from these figures. It can be observed that the TDA states [4, 0, 0], [1, 3, 0], 
[0,3, 1] (not depicted because it does not lead to any exact solution for \g\ G [0, 1]) and [0,4,0] 
do not correspond to an exact solution in the small-interaction domain, which is consistent with 
what we derived earlier. Some states are clearly connected to a given state for the whole domain 
of \g\ G [0,1] (such as e.g. [3,0,1] with the 3rd excited state), whereas other TDA states have 
no corresponding exact solution over the domain (such as e.g. [0,3, 1]). It can also occur that 
one TDA state can connect to different exact states, however for different domains of \g\ (such 
as e.g. [3, 1,0], which gives rise to the ground state for small \g\ as well as to the 2nd excited 
state for larger \g\). In addition, we can verify that the observed clustering around \g\ ~ 1 
is indeed related to the number of collective TDA eigen modes v\, consistent with the earlier 
observations. 

4. Conclusions 

In conclusion, an algebraic technique to connect the exact eigen states of the BCS Hamiltonian 
with the predictions of Tamm-Dancoff Approximation is presented. The connection is made 
via a parametrised algebra, giving rise to the quasi-spin su(2) algebra and the Heisenberg-Weyl 
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Figure 4. The TDA eigen modes and corresponding RG variables for given TDA multiphonon 
states as a function of the interaction \g\. More information is discussed in the text. 
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Figure 5. The TDA eigen modes and corresponding RG variables for given TDA multiphonon 
states as a function of the interaction \g\. More information is discussed in the text. 



hw(l) in the two limits of the parameter. Using this parametrised algebra for the Bethe Ansatz 
wavefunction, we obtain the corresponding parametrised Richardson-Gaudin equations, which 
reduce to the exact Richardson-Gaudin equations in the quasi-spin limit and the Tamm-Dancoff 
secular equation in the Heisenberg-Weyl limit. These parametrised RG equations allow us to 
make a clear-cut connection between the exact states of the BCS Hamiltonian and the TDA 
predictions. 
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